INTERIOR PENALTY DISCONTINUOUS GALERKIN METHOD ON 
VERY GENERAL POLYGONAL AND POLYHEDRAL MESHES 
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Abstract. This paper focuses on interior penalty discontinuous Galerkin methods for second 
order elliptic equations on very general polygonal or polyhedral meshes. The mesh can be composed 
of any polygons or polyhedra which satisfies certain shape regularity conditions characterized in 
a recent paper by two of the authors in I17| . Such general meshes have important application 
in computational sciences. The usual conforming finite element methods on such meshes are 
either very complicated or impossible to implement in practical computation. However, the interior 
penalty discontinuous Galerkin method provides a simple and effective alternative approach which is 
efficient and robust. This article provides a mathematical foundation for the use of interior penalty 
discontinuous Galerkin methods in general meshes. 
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1. Introduction. Most finite element metliods are constructed on triangular 
and quadrilateral meshes, or on tetrahedral, hexahedral, prismatic, and pyramidal 
meshes. To extend the idea of the finite element method into meshes employing gen- 
eral polygonal and polyhedral elements, one immediately faces the problem of choosing 
suitable discrete spaces on general polygons and polyhedrons. This issue has rarely 
been addressed in the past, partly because it can usually be circumvented by divid- 
ing the polygon or polyhedron into sub-elements using only one or two basic shapes. 
However, allowing the use of general polygonal and polyhedral elements does provide 
more flexibility, especially for complex geometries or problems with certain physical 
constraints. One of such example is the modeling of composite microstructures in 
material sciences. A well-known solution to this problem is the Voronoi cell finite ele- 
ment method [SJ [SI UHl US] , in which the mesh is composed of polygons or polyhedrons 
representing the grained microstructure of the given material. The main difficulty of 
constructing conforming finite element methods on Voronoi meshes is that, the finite 
element space has to be carefully chosen so that it is continuous along interfaces. Al- 
though the constructions on triangles, quadrilaterals, or three-dimensional simplexes 
are straight forward, it is not easy for general polygons and polyhedrons. Probably 
the only practically used solution is the rational polynomial interpolants proposed by 
Wachspress |16| . in which rational basis functions are defined using distances from 
several "nodes" . An important constraint in the construction of the Wachspress basis 
is that, the rational basis functions need to be piecewise linear along the boundary 
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of every element, in order to ensure conformity of the finite element space. This 
not only limits the approximation order of the entire Wachspress finite element space, 
but also complicates the construction. The Wachspress element has gained a renewed 
interest recently [6j[7l[TJ. However, as we have pointed out above, its construction 
is complicated and usually requires the aid of computational algebraic systems such 
as Maple. 

Another practically important issue is to define finite element methods on hybrid 
meshes. Hybrid meshes are frequently used nowadays. It can handle complicated 
geometries, and can sometimes reduce the total number of unknowns. Another pos- 
sible reason for using the hybrid mesh is that, some engineers argue that in three- 
dimensions, a hexahedral mesh yields more accurate solution than a tetrahedral mesh 
for the same geometry [181119) . as partly verified by numerical experiments. However, 
pure hexahedral meshes lack the ability of handling complicated geometries. Hence 
a hybrid mesh becomes a welcomed compromise between accuracy and flexibility. 
For conforming finite element methods based on hybrid meshes, continuity require- 
ments on interfaces must be satisfied. Such a coupling is straight-forward for the 
if ^-conforming finite elements on a triangular-quadrilateral hybrid mesh. However, 
for three-dimensional meshes, high order finite elements, or other complicated finite 
element spaces, it usually requires special treatments. 

An alternative solution, that can address both issues mentioned above, is to use 
the weak Galerkin method proposed in [17]. The weak Galerkin method uses discon- 
tinuous piecewise polynomials inside each element and on the interfaces to approxi- 
mate the variational solution. In |17j , the authors have proved optimal convergence of 
the weak Galerkin method for the mixed formulation of second order elliptic equations 
on very general polygonal and polyhedral meshes. Most of the existing error analysis 
of finite element methods assume triangular, quadrilateral, or some commonly-seen 
three-dimensional meshes. To our knowledge, it is the first time that optimal conver- 
gence for the finite element solutions has been rigorously proved in |17j for general 
meshes of arbitrary polygons and polyhedrons. 

The discontinuous Galerkin method imposes the interface continuity weakly, and 
is known to be able to handle non-conformal, hybrid meshes as well as a variety of 
basis functions. There have been many research works in this direction, for example, 
nodal discontinuous Galerkin methods (Sj El [l2] for hyperbolic conservation laws. 
However, we would like to point out that so far there has been no theoretical analysis 
on the convergence rate of discontinuous Galerkin method, on very general polygonal 
or polyhedral meshes yet. Motivated by the work in \T7\ . here we would like to 
fill the gap. The objective of this paper is to establish the theoretical analysis of 
the interior penalty discontinuous Galerkin method [2] for elliptic equations on very 
general meshes and discrete spaces. 

The paper is organized as follows. In Section 2, we briefiy describe the interior 
penalty discontinuous Galerkin method in an abstract setting. In Section 3, several 
assumptions on the discrete spaces are listed, which form a minimum requirement 
for the well-posedness and the approximation property of the discrete formulation. 
Abstract error estimations are given. In Section 4, we discuss choices of meshes and 
discrete spaces that satisfy the assumptions given in Section 3. Finally, numerical 
results are presented in Section 5. 
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2. The model problem and the interior penalty method. Consider the 
model problem 



(2.1) 




in il, 
on dn, 



where n e W^{d = 2,3) is a closed domain with Lipschitz continuous boundary, and 

For any subdomain K C SI with Lipschitz continuous boundary, we use the stan- 
dard definition of Sobolev spaces H'^{K) with s > (e.g., see [11 14] for details). The 
associated inner product, norm, and seminorms in H^{K) are denoted by (•, •)s,_fs:, 
II • \\s.K, and I • Is.if, respectively. When s = 0, H^{K) coincides with the space of 
square integrable functions L^{K). In this case, the subscript s is suppressed from 
the notation of norm, semi-norm, and inner products. Furthermore, the subscript K 
is also suppressed when K = VI. Finally, all above notations can easily be extended 
to any e C dK. For the inner product on e, we usually denote it as (•, •)£ in stead 
of (•, as it can be replaced by the duality pair when needed. 

For simplicity, we assume that £7 satisfy certain conditions such that Equation 
(|2.ip has at least H"^ regularity with r > 3/2, that is, the solution to Equation (12. ip 
satisfies u £ H^{n) and 

(2.2) \\u\\r < CM. 

This assumption is standard in the practice of interior penalty discontinuous Galerkin 
methods, as it ensures that the exact solution u also satisfies the discontinuous 
Galerkin formulation, and thus the a priori error estimation can be easily derived 
in a Lax-Milgram framework. However, such a regularity assumption is not necessary 
in the practice of interior penalty methods. A well-known technique, which was first 
proposed by Gudi [11], is to use a posteriori error estimation to derive an a priori 
error estimation for the interior penalty method, with only minimum regularity re- 
quirement u £ H^{^). We believe that the same technique applies for the general 
polygonal and polyhedral meshes, as long as a working a posteriori error estimation 
is available. However, here we choose to completely skip this issue, as it is not the 
main purpose of this paper. 

Assume that for all set K discussed in this paper, including itself, the unit 
outward normal vector n is defined almost everywhere on dK. Note this is true for 
all polygonal and polyhedral elements with Lipschitz continuous boundaries. Since 
the exact solution u G i7''(0) with r > 3/2, it is clear that for any smooth function v 
defined on K, 

(Vu, Vv)k - (Vu • n, v)qk = if, v)k, 
where (•, ■)k is the L^-inner product in L'^{K) and (•, ■)dK is the L^-inner product in 

Let 7?i be a partition of the domain i7 into non-overlapping subdomains/elements, 
each with Lipschitz continuous boundary. Here h denotes the characteristic size of the 
partition, which will be defined in details later. The interior interfaces are denoted 
by e = i?i n i?2, where Ki, K2 G Th- Boundary segments are similarly denoted 
by e = X n 9ri, where K E Th- Denote by £h the set of all interior interfaces and 
boundary segments in Th, and by 8^ = £fi\ dfl the set of all interior interfaces. For 
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every K G 7^, let \K\ be the area/volume of and for every e e S^, let |e| be its 
length/area. Denote the diameter of e £ £h and /i^^ the diameter of K ^ Th- 
Clearly, when e C dK, we have he < /i^. Finally, define h = maxKeTh to be the 
characteristic mesh size. 

Notice that Th defined above is a very general mesh/partition on f2, as we do not 
specify the shape and conformal property oi K G Th- The interior penalty discontin- 
uous Galerkin (IPDG) method can be extended to such a general mesh, without any 
modification of the formulation. However, to ensure its approximation rate, certain 
conditions must be imposed on 7/i and the discrete function spaces. In this paper, we 
are interested in discussing the minimum requirements of such conditions. First, we 
shall give the formulation of the interior penalty discontinuous Galerkin method. 

Let Vfc be a finite dimensional space of smooth functions defined on K G Th- 
Define 

Vh^{v£ L^n) : v\k e Vk, for all K e 

and 

V{h) ^Vh+I H^{n) n II H' {K) ] , where r > ^. 
\ Ken I 

For any internal interface e = K\ H K2 G S-n-, let ni and W2 be the unit outward 
normal vectors on e, associated with K\ and K2, respectively. For v G ^(/i), define 
the average {Vw} and jump \v\ on e by 

On any boundary segment e = -ftT H 917, the above definitions of average and jump 
need to be modified: 

where is the unit outward normal vector on e with respect to K . 
Define a bilinear form on V{h) x V{h) by 

A{u,v)^ (Vw,Vi;)K- ^({Vm}, H)e 

KeTh ee£h 

ee£h ee£h 

where 5 = ±1, and a > 0. when 6 = 1, the bilinear form A{-,-) is symmetric. 
The constant a is usually required to be large enough, but still independent of the 
mesh size h, in order to guarantee the well-posedness of the discontinuous Galerkin 
formulation. Details will be given later. 

It is clear that the exact solution u to Equation (|2.1|) satisfies 

(2.3) A{u, v) = (/, v) for ah veVh, 

as [u] vanishes on all e G £h- Hence the following interior penalty discontinuous 
Galerkin formulation is consistent with Equation (j2.ip : find Uh G Vh satisfying 

(2.4) A{uh,v) = {f,v) iovallveVh. 

Finally, we would like to point out that the formulation (|2.4p is computable, as 
long as each finite dimensional space Vk has a clearly defined and computable basis. 
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3. Abstract theory. Define a norm ||| ■ ||| on V{h) as following: 



(3.1) \\\v\f= ^ ||v^;||^+^/ie||{Vz;}||^ + a^^| 



Mil' 



By the Poincare inequality, ||| • ||| is obviously a well-posed norm on V{h). 

Next, wc give a set of assumptions, which form the minimum requirements guar- 
anteeing the well-posedness and the approximation properties of the interior penalty 
discontinuous Galerkin method. 

11 (The trace inequality) There exists a positive constant Ct such that for all 
K eTh &nd0 G H^{K), we have 

(3.2) \\O\\lK<CT{h],'\\0\\j, + hK\m\K)- 

12 (The inverse inequality) There exists a positive constant Cj such that for all 
K G Th, (j) &Vk and G jf-VK where i — 1, . . . ,d, we have 



(3.3) \\y<l>\\K<Cjh],^<j>\\ 



K- 



13 (The approximability) There exist positive constants s and Ca such that for 
all V € iJ''+i(fi), we have 



(3.4) inf III; - Xh\\ <Ca\Y. h'Ml+iA 



1/2 



The abstract theory of the interior penalty discontinuous Galerkin method can 
be entirely based on Assumptions 11-13. 

Lemma 3.1. Assume 11-12 hold. The bilinear form A{-,-) is bounded in V{h), 
with respect to the norm ||| • ||| . Indeed, 

1 + a 

A{u,v) < |||u||| |||w||| for all u, v £ V{h). 

Furthermore, denote Ci — Ct(1 + C/)^. Then for any constant < C < 1 and 
Q! > 4(t-c''y-'^ , the bilinear form A{-,-) is coercive on Vh- That is, 

A{v,v)>-^\\vf forallveVh. 

Proof. The boundedness of A{-, ■) follows immediately from the Schwarz inequal- 
ity. Here we only prove the coercivity. First, notice that for all w € 14, by assumptions 
11-12 and the fact that he < Hk for all e € dK n £h, 

E ^e||{Vz;}||^ < E ( E ^ellV^II^ ) 

< E f^KW^^faK 

Ken 

< E hK(cT{l + Cj)h-^'\\Vv\\j, 
Ken ^ 

Ken 
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Then, by the Schwarz inequahty, the Young's inequahty and assumptions 11-12, we 
have 



'&£h e££h ee£h 

Ken \ ee£h 
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where e is chosen to be (j^^^^^ for any given constant < C < 1. Clearly, for such 
an e, we have 1 — (1 + S)eCi = C and 

'-^-^ ^ 4(1-C)2- 
Combine the above and let a > ■^^3^^' have 



AK«)= E l|Vr^||i.-(l + <5)E({V«}, [v])e+aJ2j-\\[v]\\l 

KeTh ee£h e&£h 

>(i-(i+*>c-,) j:iiv*+(i-i±i)ox:i-iiME 

>cf E llvHI?. + «Ef IIHII^) 

^ C 2 

> \\\v\\\ . 

□ 

Lemma 13.11 guarantees the existence and uniqueness of the solution to Equation 
(|2.4p . In the rest of this paper, we shall always assume a > -^p^^c^- Let u and 
be the solution to equations (|2.ip and (|2.4p . respectively. By subtracting (I2.3P from 
(j2.ip . one gets the standard orthogonality property of the error, 

A{u — Ufi, Vfi) — for all v G Vh- 

Then clearly, for all G V/j, 

lllx/i - Uh\f < ^ A{xh - Uh, Xh - Uh) 
= — ^ — MXh -u,Xh- Uhj 

< ^ iWXh - u\\\ \\\xh - Uh\\\. 

Then, using the triangle inequality, 

|||w-u/i|||< inf I lll-u - Xftlll + lllxh - 

< 1 + r< inf ll|w-Xh|| 

\ La J xheVh 
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Combine this with assumption 13, we get the following abstract error estimation: 
Theorem 3.2. Assume 11-13 hold, C be a given constant in (0,1) and a > 

^4(i-C)'^ • ^ '^^'^ solution to equations 12.1]) and 112. 4\ ), respectively. 

Then 



\\\u-Uh\\\<CA[l + 



(l + Ci)(l + a)' ' ^ 



Ca 



(j2 h'KMl+i,K 



Finally, we derive the L? error estimation by using the standard duality argument. 
Let (5 = 1, that is, the bilinear form •) is symmetric. Consider the following 
problem 

{— A0 =- u — Uh in i7, 

(/) = on dVt. 

Here again, we assume that the domain satisfies certain condition such that has 
H"^ regularity, with r > 3/2. Let G Vh be an approximation to </) such that they 
satisfy Assumption 13. Clearly 

1|U - Uhf = (-A0, U-Uh)= J2 (^'^' " "''))^ " J2 {{'^<P}, [U - Uh])e 

= u-Uh) ^ A{(i> - 4>h, u- Uh) 

< ^—^1110 - 0/1 III III" - Uh\\\ 

a 

l + a f V^^ 

< Y1 ^/f"""^''~^'''^ll'?^llmi„{r,s+l}.K lh-"/i|l|- 



This gives the following theorem 

Theorem 3.3. Assume 11-13 hold, S = I, C be a given constant in (0,1), 

^ ^4(i~c)^ ' '^^'^ elliptic equation h2.1\) has regularity with r > 3/2. Let u 
and Uh be the solution to equations h2.1]) and {2.4^ , respectively. Then 

Wu^UhW < ii^CACfl,/i'"'"{'-i'^>|||7/-u„|||. 



4. Requirements on meshes and discrete spaces. On triangular or quadri- 
lateral meshes, the usual tool for proving assumptions 11-13 is to use a scaling argu- 
ment built on affine transformations. However, on general polygons and polyhedrons, 
it is not clear how to define such affine transformations. The assumptions 11-13 were 
first validated in [17] for general polygonal and polyhedral meshes that satisfy a set 
of conditions introduced in [17]. Such conditions can be stated as follows. 

All the elements of Th are assumed to be closed and simply connected polygons 
or polyhedrons. We make the following shape regularity assumptions for the partition 
Th. 

Al: Assume that there exist two positive constants and pe such that for every 
element K G Tk and e G f^, we have 

(4.1) p.h'^K < \K\, p,ht' < \e\. 
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A2: Assume that there exists a positive constant k such that for every element 
K eTh and e G dK Ci £h , we have 



(4.2) kHk < he- 



A3: Assume that for every K G Th, and e e dK D £h, there exists a pyramid 
P(e, iiT, yle) contained in K such that its base is identical with e, its apex is 
e K, and its height is proportional to Hk with a proportionality constant 
bounded away from a fixed positive number a* from below. In other 
words, the height of the pyramid is given by (Tf^hn such that > a* > 0. 
The pyramid is also assumed to stand up above the base e in the sense that 
the angle between the vector Xg — Ag, for any Xe G e, and the outward normal 
direction of e is strictly acute by falling into an interval [0, Oq] with 9o < Tr/2. 

A4: Assume that each K G Th has a circumscribed simplex S{K) that is shape 
regular and has a diameter hs(^ji-j proportional to the diameter of K; i.e., 
hs(K) ^ l*hK with a constant 7* independent of K. Furthermore, assume 
that each circumscribed simplex S{K) intersects with only a fixed and small 
number of such simplexes for all other elements K gTh- 

Under the above assumptions, the following results have been proved in |17j : 

Lemma 4.1. (The trace inequality). Assume A1-A3 hold on a polygonal or 
polyhedral mesh. Then II is true. 

Lemma 4.2. (The inverse Inequality) . Assume A1-A4 hold on a polygonal or 
polyhedral mesh and each Vk is the space of polynomials with degree less than or equal 
to n. Then 12 is true with Cj depending on n, hut not on hpc or \K\. 

Lemma 4.3. Assume A1-A4 hold on a polygonal or polyhedral mesh and each 
Vk is the space of polynomials with degree less than or equal to n. Let Qh be the L? 
projection onto Vh. Then for all < s < n and v G H"~^^{Q), 



Ken 
Ken 



It is not hard to see that 13 follows immediately from Lemma l473l Indeed, notice 
that as long as II and 12 are true and v e iJ'"(f2) with r > 3/2, we have 



\\\v-QHv\\\<il + C^) \\S/{v-Q,,v)rK + aY, -\\[v-Qf,v]\\l 
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-QhvfK<CQoh''^^+'^\\v\\l 
-Qhv)\\l<CQ^h'^\\v\\^,^,. 



where Ci = Ct(1 + C/)^. Next, notice that by A2, II and Lemma 




= E -j-wi^ -Q'^^)\k\\Ik 

<^(Cqo + Cqi)/z2^||z;||^+i. 

Combine the above, we have 

Lemma 4.4. (The aproximability) Assume A1-A4 hold on a polygonal or poly- 
hedral mesh and each Vk is the space of polynomials with degree less than or equal to 
n. Then for all ^ < s < n and v G there exists a constant Ca independent 

of h such that 

inf |||w-X/i||| < CAh''\\v\\s+i. 

Xh&Vh 

Here s > is added so that \\\v — Xh\\\ is well-defined. 

5. Numerical Examples. Finally, we present numerical results that support 
the theoretical analysis of this paper. We fix the coefficients S — 1 and a = 10, 
since the purpose of the numerical experiments is to examine the accuracy of the 
interior penalty discontinuous Galerkin method on arbitrary polygonal meshes, not 
for different coefficients. Consider the Poisson's equation on J7 = (0, 1) x (0, 1) with 
the exact solution u = sin{2Trx) cos(27ry). Clearly w = on dfl. For simplicity of the 
notation, we denote 

\ 1/2 

Uh\l,h = I > J N{u-Uh)\K 



The first test is performed on a non-conformal triangular-quadrilateral hybrid 
mesh. The initial mesh and the mesh after one uniform refinement are given in 
Figure 15.11 A sequence of uniform refinements are then applied to generate a set 
of nested meshes. Notice that the meshes are non-conformal and there are hanging 
nodes. However, the interior penalty discontinuous Galerkin method can deal with 
such meshes without special treatments. We solve the Poisson equation using the 
interior penalty discontinuous Galerkin formulation (12. 4p on these meshes, where the 
local discrete spaces Vk are taken to be Pi polynomials on each K G Th, matter 
whether K is a. triangle or quadrilateral. The semi-norm and the norm of the 
errors are reported in Table 15.11 and Figure 15.21 These errors are computed using a 
5th order Gaussian quadrature on triangles. For quadrilateral elements, the errors 
can be conveniently computed by dividing the quadrilateral into two triangles and 
then applying the Gaussian quadrature. Our results show that the semi-norm has 
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Fig. 5.1. Initial and refined mesh for test 1. 




Table 5.1 
Convergence rates for test 1. 



h 


1 

16 


1 

32 


1 

64 


1 

128 


1 

256 


Oih"-), r = 


\U - Uh\l,h 


1.2006 


0.5904 


0.2917 


0.1452 


0.0725 


1.0124 


\\U - Uh\\ 


0.0551 


0.0159 


0.0042 


0.0011 


0.0003 


1.9270 



an approximate order of 0{h), while the norm has an approximate order of 0{h^), 
as predicted by the theoretical analysis. 

In the second test, we consider a hybrid mesh containing mainly hexagons, but 
with a few quadrilaterals and pentagons. Indeed, it is derived by taking the dual 
mesh of a simple triangular mesh. In Figure 15.31 the initial triangular mesh and its 
dual mesh are shown. By refining the triangular mesh and computing its dual mesh, 
we get a sequence of hexagon hybrid meshes. Again, we solve the interior penalty 
discontinuous Galerkin formulation (|2.4p on these hexagon hybrid meshes, with the 
local discrete spaces Vk of Pi polynomials. The semi-norm and the norm of 
the errors are reported in Table 15.21 and Figure 15.41 Optimal convergence rates are 
achieved. 
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Fig. 5.2. Convergence rates for test 1. 
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Table 5.2 
Convergence rates for test 2. 



h 


i 

16 


i 

.'52 


i 

64 


i 

128 


i 

2. 5 6 


Oih"-), r = 




0.8139 


0.38()8 


0.1894 


0.U941 


0.047U 


1.0270 


\\U-Uh\\ 


0.0461 


0.0129 


0.0034 


0.0009 


0.0002 


1.9393 
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